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Scaling behavior is studied of several dominant eigenvalues of spectra 
of Markov matrices and the associated correlation times governing critical 
slowing down in models in the universality class of the two-dimensional Ising 
model. A scheme is developed to optimize variational approximants of pro- 
gressively rapid, independent relaxation modes. These approximants are used 
to reduce the variance of results obtained by means of an adaptation of a 
quantum Monte Carlo method to compute eigenvalues subject to errors pre- 
dominantly of statistical nature. The resulting spectra and correlation times 
are found to be universal up to a single, non-universal time scale for each 
model. 
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One of the most remarkable characteristics of critical behavior is universality. For in- 
stance, it is generally accepted that upon approach of a critical point, the correlation length 
diverges with a power law and that the exponents are universal in the sense that they depend 
only on the qualitative features characterizing the direction of approach of the critical point. 
The definition of the correlation length can be based on any of various correlation functions, 
the most obvious ones of which are the order parameter and energy autocorrelation func- 
tions. The nature of the definition shows up in the amplitude of the power law: Within a 
given universality class, this amplitude is of the form mA, where A is universal, but depends 
on the observables used in the definition of the correlation length; m is non-universal, yet 
independent of these detail, and a function only of the representative of the universality 
class. 

In terms of the spectrum of the transfer matrix one can define an infinity hierarchy of 
correlation lengths. Indeed, the spectrum of the transfer matrices, which has been studied 
extensively for two-dimensional systems, displays such universal amplitudes.0 For statics 
the transfer matrix generates translations in space. The Markov matrix is its analog for 
translations in time in for stochastic dynamics. From this perspective, the analog of the 
correlation length is the correlation time and in this Letter we address universality of the 
critical point amplitudes that describe the vanishing of gaps of the spectrum of the Markov 
matrix and the slowing-down of the independent modes. 

The dynamic critical exponent z links the divergences of the correlation length and time, 
£ and time r, via the relation r ~ £ 2 . For a finite system of linear dimension L, with 
its thermodynamic fields fixed at the critical values of the infinite system, this becomes 
t ~ L z , since the correlation length is limited by the size of the system. We denote by 
1 = Alo > Aii > • • • the eigenvalues of the Markov matrix, and by An = 1 — Xu the 
spectral gaps. For the associated correlation times one has r^ 1 = L d In A^ and the latter 
are expected to display the following scaling 
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r Ll « L~ d All -> mt A,L z , (1) 



for L — > oo. Here and z are universal, and m t is a non- universal, metric factor.13 
Furthermore, d is the dimensionality of the system, which enters because we consider a 
Markov matrix that evolves the system in a local sense only. More specifically, we study 
single-spin-flip dynamics in Ising models defined on a square lattice of linear dimension L 
with nearest and next-nearest neighbor couplings K and K' and periodic boundaries. We 
focus on models described by three ratios (3 = K'/K, namely (3 = —4,0,1, the opposite-, 
nearest-, and equivalent-neighbor models. 

The dynamics is generated by the Markov matrix P, element P(S', S) of which defines, 
given a configuration S, the conditional probability of a transition to a configuration S'. If S 
and 5" differ by more than one spin, P(S', S) = 0. If both configurations differ by precisely 
one spin, 



P(3',S) = jl-tanh 



H(S') - H(S) 
2kT 



(2) 



where Ti denotes the spin Hamiltonian. The diagonal elements P(S, S) follow from the 
conservation of probability J2s' S) = 1. 

We compute the spectrum of P by means of a method used previously for a single 
eigenstateB generalized to several dominant eigenvalues of the Markov matrix. This method 
was introduced by Ceperley and Bernu in the context of quantum Monte Carlo methods!. 
A crucial element in our approach is the construction of optimized trial states, and for this 
purpose we generalize ideas of Umrigar et aln. 

The condition of detailed balance is used to define a stochastic process that has the 
Boltzmann distribution exp[— TC(S)/kT) = i/jb(S) 2 as its stationary state. Consequently, 
P(S',S) = ^rgjr P(S', S)iPb(S) is symmetric in S and 5". We write eigenvectors of the 

transform P in the functional form ^^(S)^b(S) defined in Eqs. (12) and (13) of Ref. || 
That is, in the first place we restrict ourselves to translationally and rotationally modes, 
which are even or odd under spin inversion, and, secondly, the trial functions are written as 
linear combinations of monomials of the magnetization and other long-wavelength Fourier 
transforms of the spin configuration. 

We generalize to simultaneous optimization of multiple trial states, a powerful method0~l 
of optimizing a single many-parameter trial function. The latter is done by minimization 
of the variance of the configurational eigenvalue: Suppose that ifa(S,p) is the value of the 
trial function ?/>t for configuration S and some choice of the parameters p to be optimized. 
The configurational eigenvalue X(S,p) of a spin configuration S is defined by 

ih\S,p) = X{S,p)ih{S,p), (3) 

where the prime indicates matrix multiplication by P, i.e., f'(S) = J2s' P{^7 S')f(S') for 
arbitrary /. The optimal values of the variational parameters are obtained by minimization 
of the variance of X(S,p), estimated by means of a Monte Carlo sample. We refer to Ref. |5| 
for details and mention only one of the key features of this method: in the ideal case, i.e., 
for an exact eigenstate ip^, the variance vanishes if it were to be computed exactly but also 
if one employs an approximate Monte Carlo expression. A similar zero-variance principle 
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holds for the method of simultaneous optimization of several trial states to be discussed 
next. 

For simplicity of presentation we first generalize the above method to the more general 
ideal case in which one can exactly compute m eigenvalues of the Markov matrix P. Suppose 
we have m basis states ip^i, i = 1, ■ ■ ■ ,m and again M spin configurations S a , a = 1, . . . , M 
sampled from ^b 2 - The case we consider is ideal in the sense that we assume that these 
states ipTi span an m-dimensional invariant subspace of P. In that case, by definition there 
exists a matrix A of order m such that 

m 

i>' Ti (s a ) = Y, k ^T j (Sa), (4) 

Again, the prime on the left-hand side of this equation indicates matrix multiplication by 
P. If M is large enough, A is for all practical purposes determined uniquely by the set of 
equations (||) and one finds 

A = N~ l V (5) 

where 

Ntj = Z- 1 EtU ^Ti(S a )ip Tj (S a ) 

% = Z- l Y% = ^Ti{S a W Tj {S a ), (6) 

and where Z is an arbitrary normalization constant; again, the prime indicates matrix 
multiplication by P. In the non-ideal case, the space spanned by the m basis states ipn is 
not an invariant subspace of the matrix P. In that case, even though Eq. (f|) generically 
has no true solution, Eqs. (|5]) and (Bp still constitute a solution in the least-squares sense, 
as may be verified by solving the normal equations. 

If a set of states span an invariant subspace, so does any non-singular set of linear com- 
binations. In principle, the optimization criterion should have the same invariance. The 
spectrum of the matrix A has this property, which suggests that one subdivide the sam- 
ple in subsamples and compute the variance of the local spectrum over these subsamples. 
In practice, however, precisely this invariance gives rise to a near-singular non-linear opti- 
mization problem. Therefore, to avoid slow or no convergence, we add a contribution to 
the above least-squares merit function to ensure that the basis states themselves are good 
approximate eigenstates, rather than just their linear combinations, and we use an iterative 
optimization procedure: First a combination of the single and multi-eigenstate merit func- 
tions is used, and finally the resulting approximate eigenstates are optimized one at a time 
using the single-state procedure only. Unfortunately, this method is capricious and often we 
proceed by trial and error. 

The variational states can be used directly only to obtain results with systematic errors, 
but these can be suppressed by the quantum Monte Carlo projection method introduced by 
Ceperley and BernuQ. Define generalized matrix elements 

Nij(t) = (ifailPtlifaj) 
V l3 {t) = (V^|P t+1 |^->- (7) 
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For t = 0, Eqs. @ are Monte Carlo estimators for these matrix elements, apart from the 
inconsequential normalization constant Z. One can view the matrix elements for t > 
as having been obtained by the substitution |^ Ti ) — > P*/ 2 |^ Ti ), which implies that spectral 
weights of "undesirable" states lower down in the spectrum are reduced. The matrix elements 
in Eqs. (0) are the following time auto- and cross-correlation functions of the Markov process 
generated by the matrix P: (^Ti(5'o)'?/'Tj(5 , t))^ B 2 and (ipTi(So)ip' T j(St)) ^ , where So and St 
are spin configurations separated in time by t single-spin flips. 

It should be noted that in the limit of vanishing statistical error, each eigenvalue estimate 
obtained by the above method is bounded from above by the corresponding exact eigenvalue. 
The reader is referred to the work of Ceperley and Bernu in Ref. |6| for further details 
and references. The systematic error decreases for increasing projection time t while the 
statistical error increases. An optimal intermediate t has to be chosen, which yields biased 
estimators and some uncertainty in the reliability of statistical error estimates. 

Of the three Ising-like models investigated here, the critical point is exactly known 
only for the nearest-neighbor model, where it occurs at K — K c (0) = |ln(l + V2). 
The critical points of the two crossing-bond models — K c (l) = 0.1901926807(2) and 
Kc(— i) — 0.6972207(2) — were determined by means of a transfer-matrix technique com- 
bined with finite-size scaling^. This analysis confirmed with a high precision that the two 
crossing-bond models with belong to the static Ising universality class. 

Monte Carlo averages were taken over 1.2 x 10 8 spin configurations, for system sizes in 
the range 5 < L < 20. For the nearest-neighbor model these samples were separated by a 
number of Monte Carlo steps per spin equal to one for L = 5 and increasing quadratically 
to ten for L = 20. For the other systems these numbers where multiplied by the appropriate 
scale factors. These surprisingly short intervals are possible because the convergence of the 
eigenvalue estimates as a function of projection time t in Eqs. (|7|) is governed by lower-lying 
Markov matrix eigenvalues. These are much smaller than the largest odd eigenvalue, which 
usually determines the relaxation rate. For the system size L = 5, the Monte Carlo results 
for the largest odd eigenvalues of the three models were compared with numerically exact 
results@. The consistency of both types of results confirms the validity of our numerical 
procedures. 

As noted before for the largest odd eigenvalue of the nearest-neighbor model!, the high 
statistical accuracy of the Monte Carlo estimates of the eigenvalue is due to the accuracy of 
the approximation of the eigenvector of the Markov matrix by the optimized trial vectors. 
The present Monte Carlo results for the largest odd eigenvalues of the nearest-neighbor 
models agree with those of Ref. []. The new data are based on statistical sample smaller by 
a factor of about 7, but the current trial vectors had more variational freedom. 

For finite system sizes L we expect corrections to the leading scaling behavior tl ~ L z . 
Following Ref. ||], we assume corrections proportional to even powers of 1/L: 

r Ll ^L z J2®kiL- 2k , (8) 

k=0 

where the series is truncated at order n c . Although we cannot exclude other powers in 1/L, 
we have used Eq. ([5]) to analyze the Monte Carlo autocorrelation times. 

Results of such fits with n c = 3 are presented in Table |. The smallest systems do not 
fit Eq. (H) well for this value of n c . However, the residuals decrease rapidly when L the 
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smallest system size included in the fit is increased. The smallest acceptable value of L , as 
judged from the \ 2 criterion, is also included in Table |. 

TABLE I. Universality of the dynamic exponent z. Results of least-squares fits for the dynamic 
exponent for three Ising-like models and for five distinct relaxation modes, identified in the first 
column: ok refers to odd mode number k and ek refers to the corresponding even mode. Subsequent 
pairs of columns list Lq, the smallest system size included in the fit, and the resulting estimates 
of z((3) for three ratios f3 = K'/K. Estimated errors are shown in parentheses. To account for 
possible lack of convergence as a function of projection time t and flaws in Eq. (^), two standard 
errors are quoted. 
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The estimates of z obtained from the largest odd eigenvalues for the three models shown 
in Table | are in a good agreement mutually and also with the result z = 2.1665 (12) of 
Ref. [5| for the nearest-neighbor model. Universality of z has independently been confirmed by 
Wang and Hu0, with a level of precision in the order of 10~ 2 . The results for the largest odd 
eigenvalues are in agreement with those obtained for the other relaxation modes. Although 
the differences do occasionally amount to 3 a, we attribute this to imperfections of Eq. 
(|) and underestimation of the statistical errors of the eigenvalues themselves. Thus we 
interpret the data in Table | as a confirmation of universality of the dynamic exponent for 
different models and modes of relaxation. 

Correlation-time amplitudes were obtained from least-squares fits using Eq. (H) with z 
fixed at y , which happens to be close to the most accurate results in Table |. Subsequently, 
the non-universal metric factors m t were computed by fitting to Eq. ([I]). Defining m t (l) = 1, 
we found m t {-\) = 2.391 ± 0.002 and m t (0) = 1.5572 ± 0.0005. Table shows results of 
the fits. Figure [I] is a semi-logarithmic plot of the effective, size-dependent amplitudes 
A Li (j3) = T Li L~ z /m t derived from the spectral gaps of the Markov matrices of the opposite-, 
nearest-, and equivalent-neighbor Ising models, (5 = — ~, and 1. 

Finally, we note that if one suppresses all but the magnetization dependence of the 
optimized trial functions, the number of nodes of the resulting functions equals the number 
of the corresponding eigenvalue counted from the top of the spectrum, which is in agreement 
with the odd-even alternation shown in Tab. 0. 
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additional support from the Advanced Research Projects Agency (ARPA), the National 
Center for Research Resources at the National Institutes of Health (NIH), IBM Corporation, 
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TABLE II. Universality of correlation-time amplitudes. Results of least-squares fits for the 
finite-size amplitudes for three Ising-like models and for five distinct relaxation processes. The 
first column and the ones labeled Lq are as in Table |. The columns labeled Ai{j3) contain the 
amplitudes defined in Eq. ([[]) for three interaction ratios = K' /K with metric factors rat as 
given in the text. Estimated errors, as defined in Table ffl, are shown in parentheses. The difference 
Ai(l) — Ai(P) divided by its error is denoted by r. 
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FIG. 1. Universality of correlation-time amplitudes. Semi-logarithmic plot of the effective, 
size-dependent amplitudes Au- To separate data point for the three models ^sign/3 was added to 
all abscissae. The data collapse predicted by Eq. (jl]) was produced by fitting two metric factors, 
mt(0) and mt(— 4). Amplitudes of odd and even states alternate in magnitude. 



10 



0.1 



0.01 



•of «3f tbf. ^ ^ ^ 0^. ^ ^ ^ ^ ^ 

0?oPc^C^CC<>occ<>o<>oooc?oonoooooooo 

rVW\-^ ^ -KX -fX -(€X -(€X -(€X -t€X H€X -)€X H€X -fx H€X -fx 
Mi □^□AnAnAnAnAnAnAnAnAnAnAnAnAnA 



10 



12 



14 



16 



18 



20 



6 



and other members of the center's Corporate Research Institute. 



1 M.P. Nightingale, in Finite-Size Scaling and Simulation of Statistical Mechanical Systems, edited 
by Privman, (World Scientific, Singapore 1990), p. 287-351. 

2 M.P. Nightingale and H.W.J. Blote, J. Phys. A 16 (1983) L657; see also references therein. 
3 V. Privman and M.E. Fisher, Phys. Rev. B 30 322 (1984). 

4 H.W.J. Blote and M.P. Nightingale, Physica A 134 (1985) 274. 

5 M.P. Nightingale and H.W.J. Blote, Phys. Rev. Lett. 76 (1996) 4548. 

6 DM. Ceperley and B. Bernu, J. Chem. Phys. 89 (1988) 6316. Also see B. Bernu, DM. Ceperley, 
and W.A. Lester, Jr., J. Chem. Phys. 93 (1990) 552; W.R. Brown, W.A. Glauser, and W.A. 
Lester, Jr., J. Chem. Phys. 103 (1995) 9721. 

7 C.J. Umrigar, K.G. Wilson and JW. Wilkins, Phys. Rev. Lett. 60, 1719 (1988); C.J. Umrigar, 
K.G. Wilson and JW. Wilkins, in Computer Simulation Studies in Condensed Matter Physics, 
Recent Developments, edited by D.P. Landau K.K. Mon and H.B. Schiittler, Springer Proc. Phys. 
(Springer, Berlin, 1988). 

8 M.P. Nightingale, Computer Simulation Studies in Condensed Matter Physics, edited by D.P. 
Landau, K.K. Mon and H.B. Schiittler, Springer Proc. Phys. (Springer, Berlin, 1997). 

9 M.P. Nightingale and C.J. Umrigar, in Recent Advances in Quantum Monte Carlo Methods, 
edited by W.A. Lester, Jr. (World Scientific, Singapore, 1997). 

°M.P. Nightingale and H.W.J. Blote, submitted to Physica A (1997). 
1 F.-G. Wang and C.-K. Hu, to appear in Phys. Rev. E (1997). 



7 



